﻿/********************************************************
 *  ██████╗  ██████╗████████╗██╗
 * ██╔════╝ ██╔════╝╚══██╔══╝██║
 * ██║  ███╗██║        ██║   ██║
 * ██║   ██║██║        ██║   ██║
 * ╚██████╔╝╚██████╗   ██║   ███████╗
 *  ╚═════╝  ╚═════╝   ╚═╝   ╚══════╝
 * Geophysical Computational Tools & Library (GCTL)
 *
 * Copyright (c) 2023  Yi Zhang (yizhang-geo@zju.edu.cn)
 *
 * GCTL is distributed under a dual licensing scheme. You can redistribute 
 * it and/or modify it under the terms of the GNU Lesser General Public 
 * License as published by the Free Software Foundation, either version 2 
 * of the License, or (at your option) any later version. You should have 
 * received a copy of the GNU Lesser General Public License along with this 
 * program. If not, see <http://www.gnu.org/licenses/>.
 * 
 * If the terms and conditions of the LGPL v.2. would prevent you from using 
 * the GCTL, please consider the option to obtain a commercial license for a 
 * fee. These licenses are offered by the GCTL's original author. As a rule, 
 * licenses are provided "as-is", unlimited in time for a one time fee. Please 
 * send corresponding requests to: yizhang-geo@zju.edu.cn. Please do not forget 
 * to include some description of your company and the realm of its activities. 
 * Also add information on how to contact you by electronic and paper mail.
 ******************************************************/

#ifndef _GCTL_SPACE_FILTER_H
#define _GCTL_SPACE_FILTER_H

#include "../core.h"
#include "math_ext.h"

namespace gctl
{
    /**
     * @brief      椭圆滤波器的基函数
     */
    enum filter_base_e
    {
        Atan,
        Eps,
    };

    /**
     * @brief      二维平面上的椭圆型滤波器参数组
     */
    struct ellipse_filter_para
    {
        double mu_x; // 滤波器的中心x坐标位置
        double mu_y; // 滤波器的中心y坐标位置
        double sigma; // 短轴半径长度（以0.5倍滤波器峰值位置计算）
        double theta; // 短轴与x轴的夹角（单位为度，逆时针为正，用于旋转滤波器）
        double len_x; // x方向空间尺度规，滤波器的最终x半径长度等于sigma/len_x（旋转前）
        double len_y; // y方向空间尺度规，滤波器的最终y半径长度等于sigma/len_y（旋转前）
        double magnify; // 滤波器峰值的乘数，滤波器自身的取值范围为0 -> 1
        double rho; // 使用arctg函数构造滤波器时变化速率（一般取1），值越小变化越平缓
    };

    /**
     * @brief      二维平面上的椭圆型滤波器
     * 
     * 此函数能够计算输入坐标位置的滤波器值，调用前需要初始化其参数类型。
     *
     * @param[in]  in_x       输入x值
     * @param[in]  in_y       输入y值
     * @param[in]  elf        滤波器参数
     * @param[in]  base_func  使用arctg函数构造滤波器（默认为真），否则使用幂函数构建滤波器。
     *
     * @return     滤波器的值
     */
    double ellipse_filter(double in_x, double in_y, const ellipse_filter_para &elf, 
        filter_base_e base_func = Atan);

    /**
     * @brief 一维高斯滤波器
     * 
     * @param in_data 输入数据
     * @param out_data 输出数据
     * @param sigma 高斯分布的标准差
     * @param k_size 高斯核的大小
     */
    void gauss1d_filter(const array<double> &in_data, array<double> &out_data, double sigma = 1.5, int k_size = 3);

    /**
     * @brief 二维高斯滤波器
     * 
     * @param in_data 输入数据，大小为a1_num*a2_num
     * @param out_data 输出数据，大小为a1_num*a2_num
     * @param a1_num 输入数据第一坐标轴x的大小
     * @param a2_num 输入数据第二坐标轴y的大小
     * @param sigma 高斯分布的标准差
     * @param k_size 高斯核的大小
     */
    void gauss2d_filter(const array<double> &in_data, array<double> &out_data, int a1_num, int a2_num, 
        double sigma = 1.5, int k_size = 3);

    /**
     * @brief 三维高斯滤波器
     * 
     * @param in_data 输入数据，大小为a1_num*a2_num*a3_num
     * @param out_data 输出数据，大小为a1_num*a2_num*a3_num
     * @param a1_num 输入数据第一坐标轴x的大小
     * @param a2_num 输入数据第二坐标轴y的大小
     * @param a3_num 输入数据第三坐标轴z的大小
     * @param sigma 高斯分布的标准差
     * @param k_size 高斯核的大小
     */
    void gauss3d_filter(const array<double> &in_data, array<double> &out_data, int a1_num, int a2_num, int a3_num, 
        double sigma = 0.5, int k_size = 3);

    /**
     * @brief 一维数组各个计算点对应的滑动窗口内的数据索引
     * 
     * @param data_len 一维数据总个数
     * @param win_len 窗口大小
     * @param cur_data_idx 当前计算点的索引
     * @param win_data_idx 计算对应的窗口内数据索引
     */
    void moving_window_1d(int data_len, int win_len, int cur_data_idx, array<int> &win_data_idx);

    /**
     * @brief 规则网格滑动平均
     * 
     * @param in_data 输入数据
     * @param out_data 输出数据
     * @param k_size 窗口的大小
     */
    void moving_average_1d(const array<double> &in_data, array<double> &out_data, int win_size = 3);
    
    /**
     * @brief 规则网格滑动平均
     * 
     * @param in_data 输入数据，大小为a1_num*a2_num
     * @param out_data 输出数据，大小为a1_num*a2_num
     * @param a1_num 输入数据第一坐标轴x的大小
     * @param a2_num 输入数据第二坐标轴y的大小
     * @param k_size 窗口的大小
     */
    void moving_average_2d(const array<double> &in_data, array<double> &out_data, int a1_num, int a2_num, int win_size = 3);

    /**
     * @brief 规则网格滑动平均
     * 
     * @param in_data 输入数据，大小为a1_num*a2_num*a3_num
     * @param out_data 输出数据，大小为a1_num*a2_num*a3_num
     * @param a1_num 输入数据第一坐标轴x的大小
     * @param a2_num 输入数据第二坐标轴y的大小
     * @param a3_num 输入数据第三坐标轴z的大小
     * @param win_size 窗口的大小
     */
    void moving_average_3d(const array<double> &in_data, array<double> &out_data, int a1_num, int a2_num, int a3_num, int win_size = 3);

    /**
     * @brief 使用趋势分析法计算输入网格的趋势网格
     * 
     * @param data 输入的二维矩阵
     * @param m_RowNum 滑动窗口的行数
     * @param m_ClmNum 滑动窗口的列数
     * @param m_OrderNumRow 滑动窗口的行的多项式阶次
     * @param m_OrderNumClm 滑动窗口的列的多项式阶次
     */
    void trend_2d(matrix<double> &data, int m_RowNum, int m_ClmNum, int m_OrderNumRow, int m_OrderNumClm);
}

#endif // _GCTL_SPACE_FILTER_H